
#include <iostream>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;

int main () 
{

    double dt = 0.0001;
//
// Middle line construction
//
    int Npoints = 10;
    floatTinyArray_t line[Npoints];
    double left_extreme = 0.1;
    double right_extreme = 0.9;

    double factor = (right_extreme - left_extreme) / (Npoints - 1);
    
    cout << "\n ----- Line construction ----- \n\n";

    for(int i = 0; i < Npoints; i++) {
	line[i](0) = left_extreme + i * factor;
	line[i](1) = 0.5;
	cout << " line[" << i << "] = " << line[i] << endl;
    }

    cout << "\n ----- Tracking: begins ----- \n\n";

    StructuredMesh< Uniform<double, 2> > mesh(1, 256, 1, 256);
    floatTinyArray_t deltas = mesh.getDeltas();

// Velocity field at time t 
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());

// Velocity field at time t + dt
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());

    ParticleTracking< RungeKutta4< Linear <double, 2> > > particle;

// Read the velocity at time t from a file
    InOut::readFromFile_1(u, v, 151, "./Velocity/velc.");

// Read the velocity at time t + dt from a file
    InOut::readFromFile_1(un, vn, 152, "./Velocity/velc.");

// Now track each particle on the line
    for(int i = 0; i < Npoints; i++) {
	line[i] = particle.track(line[i], deltas, u, v, un, vn, dt);
	cout << " line[" << i << "] = " << line[i] << endl;
    }

    cout << "\n ----- Tracking: ends ----- \n\n";

    return 0;
}



